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Abstract 


Introduction 


A time accurate approximate factorization 
(AF) algorithm is formulated for solution of the 
three dimensional unsteady transonic small - 
disturbance equation. The AF algorithm consists 
of a time linearization procedure coupled with a 
Newton iteration technique. Superior stability 
characteristics of the new algorithm are 
demonstrated through applications to steady and 
oscillatory flows at subsonic and supersonic 
freestream conditions for an F-5 fighter wing. 
For steady flow calculations, the size of the 
time step is cycled to achieve rapid conver- 
gence. For unsteady flow calculations, the AF 
algorithm is sufficiently robust to allow the 
step size to be selected based on accuracy 
rather than on stability considerations. 
Therefore, accurate solutions are obtained in 
only several hundred time steps yielding a 
significant computational cost savings when 
compared to alternative methods. 


Nomenclature 

c airfoil chord 

c r wing reference chord 

C p pressure coefficient 

★ 

Cp critical pressure coefficient 

Cp unsteady pressure coefficient 

normalized by oscillation amplitude 
f function defining position of wing 

k reduced frequency, wc r /2U 

M freestream Mach number 

t time, nondimensional i zed by freestream 

velocity and wing reference chord 
U freestream velocity 

oq mean angle of attack 

amplitude of pitch oscillation 
Y ratio of specific heats 

r circulation 

At nondimensional time step 

n fractional semispan 

A leading edge sweep angle 

<P disturbance velocity potential 

w angular frequency 

Subscripts 

j index of grid points in spanwise 

direction 

J index of grid point at wing root 

k index of grid points in vertical 

di recti on 

TE trailing edge 


Presently, considerable research is being 
conducted to develop finite-difference computer 
codes for calculating transonic unsteady 
aerodynamics for aeroelastic applications. 1 The 
computer codes are being developed to provide 
accurate methods of calculating unsteady 
airloads for the prediction of aeroelastic 
phenomena such as flutter and divergence. The 
most fully developed U.S. code for transonic 
aeroelastic analysis of isolated planar wings is 
XTRAN3S. 2 The code uses an alternating- 
direction implicit (ADI) finite-difference 
algorithm to calculate steady and unsteady 
transonic flows in a time-marching fashion. 
Several terms of the ADI algorithm are treated 
explicitly, though, which leads to a time step 
restriction based on numerical stability 
considerations. Experience with the code has 
shown that for applications to practical wings 
with moderate to high sweep and taper, very 
small time steps are required for the algorithm 
to be numerically stable. 3 * 6 This stability 
restriction typically results in thousands of 
time steps required to obtain converged 
solutions, which generally is many more time 
steps than is necessary to resolve the physics 
of the problem. Such solutions are 

computationally expensive, and thus, aeroelastic 
applications of XTRAN3S have generally been 
limited. An algorithm Is therefore desired 
which is robust in comparison with the ADI 
algorithm to allow for efficient application to 
aeroelastic problems. 

The purpose of this paper is to describe 
the development of a time-accurate approximate 
factorization (AF) algorithm for solution of the 
unsteady transonic smal 1 -di sturbance (TSD) 
equation. The AF algorithm involves a local 
time linearization procedure coupled with a 
Newton iteration technique which is based on 
work recently reported by Shankar, Ide, 
Gorski, and Osher 7 and Shankar and Ide. 8 In 
Refs. 7 and 8, these concepts were developed for 
the ful 1 -potential equation and the resulting 
algorithm was shown to be very robust for 
application to either steady or oscillatory 
transonic flow problems. The objectives of the 
present research were to: (1) develop a similar 

AF algorithm for solution of the unsteady TSD 
equation, (2) validate the new algorithm by 
making comparisons with available experimental 
data as well as with results computed using the 
XTRAN3S ADI algorithm, and (3) investigate the 
robustness and efficiency of the new solution 
procedure. The paper presents a description of 
the AF algorithm along with results and 
comparisons which assess this new capability. 
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Governing Equations 

In this section the equations governing the 
flow are briefly described including the ISO 
equation and boundary conditions. Oetalls of 
the coordinate transformation to map the 
physical domain into the computational domain 
also are given. 


Trailing wake: [* ] = O' (5f) 

[$x + * t J - 0 ( 5g ) 


where [ ] indicates the jump in the indicated 
quantity across the wake. The wing flow- 
tangency boundary condition is 


TSP Equation 

The flow is governed by the modified TSD 
equation which may be written in conservation 
law form as 


*2 ‘ f x‘ * 't < 6 > 

which is imposed at the mean plane of the wing. 
The plus and minus superscripts indicate the 
upper and lower wing surfaces, respectively. 


3f 


0 


+ 



at 3x ay az 


0 


( 1 ) 


Coordinate Transformation 


where 

f Q = - A» t - (2a) 

^1 “ E$ x + + G^ 2 (2b) 


The finite-difference grids in both the 
physical and computational domains are contained 
within rectangular regions and conform to the 
wing planform. Regions in the physical domain 
such as the swept and/or tapered wing are mapped 
into rectangular regions in the computational 
domain using the shearing transformation 


" * y + H Vy ( 2c ) 

f 3 - * 2 (2d) 


The coefficients A, B, and E are defined as 


A = M 2 , B = 2M 2 , E = 1 - M 2 


(3) 


Several choices are available for the 
coefficients F, G, and H depending upon the 
assumptions used in deriving the TSD equation. 2 
In this study, the coefficients are defined as 


5 = S(x,y), n = y, C * z (7) 


where 5, n, and c are the nondimensional 
computational coordinates in the x, y, and z 
directions, respectively. The TSD equation 
(Eq. (1)) may then be expressed in computational 
coordinates as 



F = - 2 (y * 1)M 2 (4a) 

G =* \ (Y - 3)M 2 (4b) 

H = - (y - 1)M 2 (4c) 


Boundary Conditions 


Boundary conditions 
field are 

imposed upon 

the flow 

Far upstream: 

<p =0 

(5a) 

Far downstream: 

<?> + = 0 
x t 

(5b) 

Far above and below: 

*z =° 

(5c) 

Far spanwise: 

X* 

II 

o 

(Sd) 

Symmetry plane: 

o 

II 

(Se) 



Alternating-Direction Implicit Algorithm 

An alternating-di rection implicit algorithm 
was developed in Ref. 2 for solution of the 
modified TSD equation (Eq. (I)). This algorithm 
forms the basis of the XTRAN3S computer code 
which can be used to calculate steady and 
unsteady transonic flow fields about planar 
wings including aeroelastic deformation 
effects. The program is capable of treating 
either forced harmonic or aeroelastic transient 
type motions. Details of the algorithm 
development and solution procedure are given in 
Ref. 2. The original XTRAN3S code used a simple 
global shearing transformati on which restricted 
applications to wings with mild sweep and taper. 
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Alternative coordinate transformations 3 ~ 5 were 
subsequently developed to allow application to 
practical wings. Although the code now can be 
applied to more general planforms, several of 
the terms of the TSO equation are treated 
explicitly in the ADI algorithm which leads to a 
time step restriction based on numerical 
stability considerations. For- application to 
wings with moderate to high sweep and taper, 
very small time steps are required for the 
algorithm to be numerically stable. For 
example, a summary of selected time steps and 
number of steps per cycle is listed in Table 1, 
the values of which are taken from Refs. 3-6. 
Cases are tabulated for three wings which 
include the F-5 wing, 3 » 4 the RAE tailplane, 5 
and the B-I wing. 6 (The calculations for these 
wings were performed using different 
computational grids and different coordinate 
transformations. ) As shown in the table, the 
stability restriction typically results in 
thousands of time steps required to obtain 
converged steady-state solutions and thousands 
of steps per cycle of forced harmonic or 
aeroelastic motion. These solutions are 
computationally expensive, and thus, aeroelastic 
applications of XTRAN3S have generally been 
limited. 

To investigate the stability problem of the 
ADI algorithm, a numerical stability analysis 
was performed for the F-5 wing. Shown in Fig. 1 
are the resulting stability boundaries 
calculated using the algorithm started from an 
initial undisturbed flow. The boundaries were 
determined by first selecting a step size of At 
= 0.05 and running the algorithm until the 
solution diverged causing program failure. The 
step size was then successively reduced and the 
calculations repeated until the boundary was 
determined. The analysis was performed for 
freestream Mach numbers of M = 0.6, 0.8, 0.9, 
and 0.95. As shown in the figure, the ADI 
algorithm is unstable in the region above and to 
the right of the boundaries; the algorithm is 
stable in the region below and to the left. The 
stability boundaries indicate that a very small 
step size is required to obtain stable solutions 


using the ADI algorithm. In fact, this step 
size is excessively small as demonstrated by the 
following example. The number of steps per 
cycle N, for oscillatory flow, is determined for 
a given time step by 


If the maximum allowable time step is 0.01 and 
the reduced frequency k is 0.1, then Eq. (9) 
indicates that over 3000 steps per cycle are 
-required for numerical stability. This is an 
order of magnitude more steps than is necessary 
for accurate aeroelastic calculations. 

To further investigate the problem, 
stability boundaries were obtained for several 
leading edge sweep angles by shearing the F-5 
planform aft. Values selected for the sweep 
angle include A = 32°, the actual sweep of the 
F-5 wing, and A - 45° and 60°. Shown in Fig. 2 
are the resulting stability boundaries for the 
freestream Mach number of 0.9. The boundaries 
indicate that as the sweep increases, 
smaller time steps are required for numerical 
stability. For the A = 60° case, for example, 
the maximum allowable time step is approximately 
0 . 002 . 


Approximate Factorization Algorithm 

An approximate factorization algorithm was 
developed as an alternative solution of the 
modified TSD equation (Eq. (I)). In this 
section, the AF algorithm is described in detail 
including the mathematical formulation, 
finite-difference discretization, boundary 
conditions, and solution procedure. 


General Description 

The AF algorithm consists of a time 
linearization procedure coupled with a Newton 


Table 1 Summary of selected time steps and steps per cycle using 
the XTRAN3S ADI algorithm reported in Refs. 3-6. 


Reference 

Wi ng 

Leading 

edge 

sweep 

angle 

Taper 

ratio 

Aspect 

ratio 

Steady 

Unsteady 

At 

total 

steps 

case 

steps/cycle 

3 

F-5 

31.9° 

0.28 

3.16 

0.01 

2000 

40 Hz, pitch 

1200 

4 

F-5 

31.9° 

0.28 

3.16 

0.01 

4000 

20 Hz, pitch 

2160 

5 

1 

RAE 

50.2° * 

0.27 

2.41 

0.0075 



1500 

70 Hz, pitch 
33 Hz, pitch 

1000 

2000 

— 

6 

B-l 

67.5° 

0.38 

1.85 

0.002 

not 

reported 

aeroelastic 

6000 
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★ = 0 


( 11 ) 


Iteration 

at 

failure 



Fig. 1 Numerical stability boundaries of the 
ADI algorithm computed from an initial 
undisturbed flow past the F -5 wing. 


Iteration 

at 

failure 



Time step. At 


Fig. 2 Numerical stability boundaries of the 
ADI algorithm computed from an initial 
undisturbed flow past the F -5 wing for 
various sweep angles. 


iteration technique. For unsteady flow 
calculations, the solution procedure involves 
two steps. First, a time linearization step 
(described below) is performed to determine an 
estimate of the potential field. Second, 
Newton Iterations are performed to provide time 
accuracy. More specifically, the TSD equation 
(Eq. ( 1 )) is written in general form as 


R(* n+1 ) = 0 


( 10 ) 


where 9 represents the unknown potentials at 
time level (n+l). The Newton iteration solution 
to Eq. ( 10 ) is then given by the first order 
Taylor series 


R(<f*) 


♦ <$> 

f= 4 > 


* 

In Eq. ( 11 ), 9 is the currently available value 
n+l n+l * 

of 9 and A<t> = - 9 , During convergence 

of the Newton iteration procedure, a? will 
approach zero so that the solution will be given 
n+l + 

by = $ . In general, only one or two 

iterations are required to achieve acceptable 
convergence since the Newton iteration process 
is quadratically convergent. For steady flow 
calculations, Newton iterations are not used 
since time accuracy is not necessary when 
marching to steady -state. 


Mathematical Formulation 


The AF algorithm is formulated by first 
approximating the time derivative terms ($tt 
and 9 xt terms) by second-order accurate 
finite-difference formulae. By substituting 
★ 

9 = 9 + A9 into the TSD equation and neglecting 

squares of derivatives of A9 (which is 
equivalent to applying Eq. ( 11 ) term by term), 
each term of Eq. ( 1 ) may then be rewritten as 


3f o 

at 


0 * r n , n-1 n-2 

A 2<j> - 69 + +9 - 9 

- A £ 

Al 


- V* 

Al 


- B 


- * . n n-1 

3 <f> - 49 + 9 V 


2 At 


3B A 

2 af 


(12a) 


!li , * (E ,* + f /2 + 6 ,*2 j 

3 x 3 x v y x Y x y ' 


+ 37 ( Ba ^x + 2F<p x A *x + 2G<f y A? y) 


(12b) 


^2 a * * * 

3y 3y ( $ y + H *x *y ^ 


+ (A9 + H9 A9 + H9 A9 ) 

3y V Ty Y x Yy Yy y x ) 

sr - It 


(12c) 

(12d) 


Summing these four terms and rearranging gives 


2 A 

At 


3 B 

2 + 2 At A *x 


ax 


(EA9 + 2 F 9 A 9 + 2 G 9 u A9 ) 


3 ★ * 

— (A9 + H9 A9 + H9 A9 ) 

3y \ 'y y x T y y y *x' 


- — ( A9 ) = 

az v V 


% 


V 
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> 


♦ 


t 


- A 


2±_ - 5» n + 44> n -I 
At^ 


i.n-2 


R = -5, 


At 2 r A 2$* - 5$ n + A*"’ 1 
x ZK 1 " £ 


Ji-2 


At 


- B 


34* - 4$ n + 4> n * 1 
Y x Y x Y x 

TEt 


D 2 \ : 4 *" + £! 
0 2At 


a * *? *? 

+ 37 ( F * x ♦ F* x + G* y ) 

a * * ★ 

+ 3j (♦y + >V y ) 

♦ If (<) 03) 


The right-hand side of Eq. (13) is simply the 
TSD equation which may be evaluated using 

known potentials <t> , $ n , , and <p n “^ . The 

left-hand side of Eq, (13) consists of terms 
containing A<j> and its spatial derivatives. At 
* n+1 

convergence $ approaches <p so that A$ 

vanishes and the left-hand side becomes zero. 

Equation (13) is transformed into 
computational coordinates using Eq. (7) and is 
rewritten in conservation form for solution by 
approximate factori zation. Each term is 
multiplied by At 2 /(2A) and the left-hand side 
is approximately factored into a triple product 
of operators yielding 



Equation (14) is solved using three sweeps 
through the grid by sequentially applying the 


operators L^, L n , 

and L c as 


C - sweep: 

A$ = - R 

(16a) 

n - sweep: 

A$ * A$ 

(16b) 

C - sweep: 

A$ = A$ 

(16c) 


L^L n L^A<P= - R(<J> > 4> n * 4» n ^) (14) 


where 


3B 


3_ , At 2 _3_ p _3_ 

35 ■ S IK 35 1 35 


l + !5f5 x At-ih--5„'Si- — Fi — 


1,1 , At 2 3 F 3_ 

L n 1 • 5 x 2* TFT F 2 1* 


i 1 C At 2 3 c 3 

L C = 1 ' 5 x JK 3; F 3 3c 


(15a) 

(15b) 

(15c) 


Fj * E5 X + 2F5 x 4> 5 + 2G5 y { £ y + 4 n ) 


Spatial Discretization 

Central difference formulae are employed 
for all of the derivatives on the left-hand 
sides of Eqs. (16) except for the second term in 
the operator (from the <*> x t term) which is 
backward differenced to maintain stability and 
the third term in the operator which is 
split into streamwise and spanwise components. 
The resulting terms are centrally differenced at 
subsonic points and the streamwise terms are 
upwind biased at supersonic points using the 
Murman 9 type-dependent mixed difference 
operator. The terms on the right-hand side of 
the 5 - sweep are also approximated using 

central -difference operators except for the 
$xt term wh ich is backward differenced and for 
terms in the streamwise direction which are 
upwind biased at supersonic points. For 
example, the term 



I 5 [ E S(*5 + H l\ 2 * G( V5 + O 2 ) (17a) 

is rewritten as 

°5 ( E5 x* 5 + ^x^ + G S ( V5 + V 2 1 

+ S 5 ^ G N ^ ^ 5 + * 2 i ( 17b ) 


5 


where G$ and Gn are the streamwise and 
spanwise components of the constant G, 
respectively, defined as 


G s = 1 - M 2 , G n = | (y - 1) M 2 - 1 (18) 


In Eq. (17b), 6 ^ is a central difference 

operator and is the Murman 9 type-dependent 

mixed difference operator that depends upon the 

sign of E + 2F$ X . 


the boundary conditions. In the following 
subsections, the numerical treatment of each of 
the boundary conditions is described in detail. 


Wing . - The wing f 1 ow-tangency boundary 
condition (Eq. ( 6 )) is imposed within the 

differencing of the term which appears 

in R (Eq. (15g)) on the right hand-side of Eq. 
(14) and within the L c operator in the time- 
linearization step for unsteady calculations. 

In general, the 3(i> )/3c term is approximated by 


In the present coding of the AF algorithm, 
the time-derivatives are implemented for 
variable time stepping to allow for step-size 
cycling to accelerate convergence to 
steady-state. In these calculations the step 
size is cycled using a standard geometric 
sequence. Also, since the l^, L n , and L c 
operators only contain derivatives in their 
respective coordinate directions, all three 
sweeps may be vectorized. This is in contrast 
with the ADI algorithm where only the streamwise 
sweep is vectori zable. 


Time-Li neari zation Step 

An initial estimate of the potentials at 
time level (n+1) is required to start the Newton 
iteration process. This estimate is provided by 
performing a time-linearization calculation. 
The equations governing the time-1 ineari zati on 
step are derived in a similar fashion as the 
equations for Newton iteration. The only 
difference is that the equations are formulated 
by linearizing about time level (n) rather than 
the iterate level (*). So by substituting 

<J> * <P n + into the TSD equation (Eq. (1)) and 
neglecting squares of derivatives of Aq>, the 
time -1 i neari zation step may be written as 


4 4 4 s o' 1 ' 1 - » n ' 2 ) (19) 


where the operators l ? , L n , and lr are 
defined by Eqs. (15a), (15b), and (15c), 

* n 

respectively, with ? replaced by <j> in 

the definitions of Fj and F 2 (Eqs. (15d) and 
15(e)). 


Boundary Conditions 

The boundary conditions are numerically 
imposed by redefining the L 5 , L n , and lr 
operators in Eq. (14) as well as the right-hand 
side R, at the appropriate grid points. The 
equation to be solved at boundary grid points 
may then be written symbolically as 


4 4 4 * - R 


( 20 ) 


where the "tilde” indicates that the quantity 
has been modified or rewritten to account for 


iJOk =- 

H C c 


k +1 


*k-l 


s k+l/2 


0 * ) ( 21 ) 
Sc -1/2 


where the derivatives on the right-hand side are 
written about half-node points. The wing is 
located equi di stant ly between grid lines so that 

in the plane directly above the wing the 

★ 

^k- 1/2 der1vative EQ* ( 21 ) is replaced 

by (f* + f t ) n+ * and ™ the P^ ane below the wing 
the ^ + 1/2 derivative in Eq. (21) is 

replaced by (f~ + f t ) n + *. The L ^ operator is 
modified in the time-linearization step since A 4 > 
is defined then to be equal to $ 0 + 1 - <*> n . 

Analogous to the 3(<}>^ )/K modifications, the 
A *^k-l /2 term replaced by 

(f* + f t ) n + 1 * + \ at 9 nd P° ints the 

plane above the wing, and the ^^ 1/2 term 

in L c is replaced by (f“ + f t ) n+1 - ( f ” + 

at grid points in the plane below the wing. 
Since these terms are known quantities, they are 
brought to the right-hand side of Eq. (16c) 
which results in a bidiagonal c - sweep 
equation. 

Wake. - The wake boundary condition (Eqs. 
( 5 f) and ( 5 g ) ) is similar to the wing 
f low-tangency boundary condition in that it is 

imposed within the 3(<p^)/ac term and the 

right-hand side of the c - sweep. The wake 
circulation r is calculated from Eq. (5g) which 
is equivalent to 


4 + 4 s 0 


( 22 ) 


Starting from the trailing edge value given by 

r JE = - <pj£, the circulation is convected 

downstream by integrating Eq. (22) using 

second-order accurate finite-difference 

approximati ons for r x and r^. The wake 
circulation is incorporated within the solution 
procedure by requiring that the perturbation 
velocity in the vertical direction be continuous 
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* 


across the wake (Eq. (5f)). This condition is 
satisfied by defining 


* * *, 
f , \ ' ( \-l * r } 

c k-l/2 S " c k-l 


(23a) 


at grid points in the plane above the wake, and 
by defining 


★ * • 
f = ^k+1 ~ r ^ ~ *k 
; k+l/2 ; k+l " c k 


(23b) 


at grid points in the plane below the wake. 
Similar substitutions are made for terms in the 

★ 

L operator by replacing $ in Eqs. (23) by A$ 

£ 

* n + 1 * 

and by replacing r in Eqs. (23) by r - r . 
Since the circulation terms are known 
quantities, they are brought to the right-hand 
side of Eq. (16c) which results in a modified 
C - sweep equation. 


Symmetry plane. - The symmetry condition 
(Eq. ( 5e ) ) is imposed at the plane of the wing 
root j * J by requiring that 


+ Vj-1/2 * ‘ ^ + ^ J +1/2 (24a * 


and 


(yj + Oj = 0 {24b) 

Equations (24a) and (24b) affect terms which 
appear on the right-hand side of Eq. (14). 
Equation (24b) causes two terms in (Eq. 
(15d) ) to vanish. The l n operator (Eq. (15b)) 
is also modified at j * J by requiring that 


(F 2 3n A ^J-l/2 * ' (F 2 3n A ^J+l/2 


(25) 


which results In an upper bidiagonal n - sweep 
equation at the wing root. The condition at the 
far spanwise boundary (Eq. (5d)) is identical to 
the symmetry condition (Eq. (5e)) and thus is 
treated similarly. 


Farf ield. - The farfield boundary 
conditions (Eqs. (5a), (5b), and (5c)) are 
imposed by writing finite-di fference 
approximations for these equations, casting them 
in the form of Eq. (20), and including them with 
the system of simultaneous equations which 
results from Eq. (14) for solution. 


Results and Discussion 


Calculations were performed for the F-5 
wing to assess the accuracy and efficiency of 
the AF algorithm. The wing has a full -span 
aspect ratio of 3.16, a leading edge sweep angle 
of 31.9°, and a taper ratio of 0.28. The 
airfoil section of the F-5 wing is a modified 
NACA 65A004.8 airfoil which has a drooped nose 
and is symmetric aft of 40% chord. The AF 
results are compared with parallel calculations 
performed using the ADI algorithm. In these 
calculations identical grids and coordinate 
transformations (Eq. (7)) were used to allow for 
a direct comparison between AF and ADI results. 
The grid and transformati on were identical to 
those used in Ref. 4. Furthermore, the results 
are compared with the experimental data from an 
F-5 wing model tested by Tijdeman, et al. 10 In 
Ref. 10, steady and oscillatory pressure 

distributions were measured at eight span 
stations along the wing. The stations were 

located at n = 0.18, 0.36, 0.51, 0.64, 0.72, 
0.82, 0,88, and 0.98. In this paper, 

comparisons are presented at the first, third, 
fifth, and seventh span stations. The tests 
covered the Mach number range from M = 0.6 to 
1.35. In this study, two Mach number cases were 
selected for application of the new algorithm. 
The first case, herein referred to as Case 1, 
was chosen to have the same freestream 

conditions as investigated in Refs. 3 and 4. In 
Case 1, the freestream Mach number was 0.9 and 
the mean angle of attack was 0°. At these 
conditions, both steady and unsteady results 
were obtained. The unsteady calculations were 
performed for the rigid wing pitching 

harmonically about a line perpendicular to the 
root at the root midchord. The second case, 

herein referred to as Case 2, was chosen to 
assess the performance of the AF algorithm for 
supersonic freestream conditions. In Case 2, 
the freestream Mach number was 1.1 and the mean 
angle of attack was 0°. No ADI calculations 
were attempted for this case since the algorithm 
is unstable for supersonic freestream cases for 
swept wings. 


Case 1: M = 0.9 and a = 0° 

o 


For Case 1, steady calculations were 
performed using both the ADI and AF algorithms. 
The ADI results were obtained using a constant 
step size of At = 0.01 which is the same as that 
reported in Refs. 3 and 4. These calculations 
were performed for a total of 4000 time steps. 
The AF results were obtained by cycling the step 
size through a range of values between At = 0.05 
and 5.0. A total of 400 time steps were run. A 
comparison of steady -state convergence between 
the ADI and AF algorithms is shown in Fig. 3. 
After 4000 time steps, the maximum |a<J)| in the 
ADI solution has been reduced to approx i mately 
10-6.5 whereas the AF solution has achieved 
similar convergence in less than 400 steps. The 
total lift and moment in the AF solution were 
converged after approximately 100 steps. Steady 
calculations were also performed for the high 
sweep case of A - 60°, to further assess the 
convergence characteri sti cs of the two 
algorithms. The ADI results were obtained using 
a constant step size of At = 0.002 and the AF 
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log 

( max|A®| ) 



Fig. 3 Comparison of steady -state convergence 
between the ADI and AF algorithms for 
the F-5 wing at M = 0.9 and Oq * 0°. 


log 

I max|AO| ) 



Time steps 


Fig. 4 Comparison of steady-state convergence 
between the ADI and AF algorithms for 
the F-5 wing sheared to A = 60° at 
M = 0.9 and Qq = 0°. 


results were obtained by cycling the step size 
through a range of values between At * 0.05 and 
0.5. A comparison of the steady-state 
convergence histories between the two algorithms 
is shown in Fig. 4. The ADI solution converges 
very slowly such that after 4000 steps, the 
maximum |a<j>I has been reduced to only 

approximately 10" 5 . The AF solution, however, 
has converged very quickly to maximum |a$| < 
10-7 i n on iy 250 steps. ' 1 

A comparison of the AOI and AF steady 

pressure distributions for the F-5 wing is shown 
in Fig. 5. These pressure distributions 
indicate that there is an embedded supersonic 

region on the wing upper surface and the flow is 
slightly supercri tical along the lower surface. 
The AF pressures are nearly identical to the ADI 
pressures and both sets of results are in good 



x/c 


Fig. 5 Comparison of F-5 wing steady pressure 
distributions at M = 0.9 and oq = 0°. 


general agreement with the experimental data. 
The AF algorithm thus produced the same 
steady-state solution for this case as the ADI 
algorithm, at one-tenth of the computational 
cost. 


Unsteady results were obtained using both 
algorithms for the wing oscillating with 
amplitude aj = 0.109° at a reduced frequency 
of k 0.137, for comparison with the 

experimental data. The ADI calculations were 
performed using 2000 steps per cycle of motion 
which corresponds to a step size of At = 0.0115, 
(In Ref. 4, similar calculations were performed 
using 2160 steps per cycle with a step size of 
At = 0.0106.) For unsteady calculations with 
the ADI algorithm, the step size is typically 
selected based on numerical stability 

considerations rather than on accuracy. 8ecause 
of the stability restriction inherent in the ADI 
algorithm, this generally leads to using many 
more time steps per cycle than is necessary to 
resolve the physics of the problem. The AF 
algorithm, however, allows the step size to be 
selected based on accuracy rather than on 
stability. Consequently, a convergence study 
was performed using the AF algorithm to 
determine the largest step size (fewest number 
of steps per cycle) that would produce converged 
answers. Unsteady results were obtained for 
100, 200, 300, and 400 steps per cycle which 
required At = 0.2293, 0.1147, 0.0764, and 

0.0573, respectively. Only one Newton iteration 
per time step was required to satisfy the 
convergence criteria of maximum I Ag>J < 10~ 6 
for the latter three cases. For thfe 100 steps 
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per cycle case, two Newton iterations per time 
step were required during approximately 60% of 
the cycle of motion. Selected results from this 
convergence study are plotted in Fig. 6 for 200 
and 300 steps per cycle. The unsteady pressure 
coefficients along the upper surface are shown 
in Fig. 6(a); the unsteady pressure coefficients 
along the lower surface are shown in Fig. 6(b). 
These coefficients are plotted as real and 
imaginary components corresponding to the 

in-phase and out-of-phase unsteady pressure 
distributions normalized by the amplitude of 
motion. The calculation for 100 steps per cycle 
of motion produced reasonable results but fairly 
large differences were observed in comparison 
with the 200 steps per cycle calculation. As 
shown in Fig. 6, the results for 200 and 300 
steps per cycle are very similar with the 
largest differences occurring near the leading 
edge and in the shock pulse region on the upper 
surface. The pressure distributions computed 
using 400 steps per cycle are identical, to 
plotting accuracy, to those computed using 300 
steps per cycle and therefore are not shown. A 
converged solution hence requires approximately 
300 steps per cycle of motion for this case, 
although the results computed using 200 steps 
may be acceptable for engineering purposes. 

Figure 7 shows a comparison between the AOI 
and AF unsteady pressures for Case 1 along with 
the experimental data. Upper surface pressure 
distributions are shown in Fig, 7(a); lower 
surface pressure distributions are shown in 
Fig. 7(b). The two sets of calculated results 
are very similar, and both sets generally agree 
well with the experimental data. The AF 
algorithm (using 300 steps per cycle and one 
Newton iteration per time step) thus produced 
approximately the same unsteady solution for 
this case as the AOI algorithm (using 2000 steps 
per cycle), at 30% of the computational cost. 


Case 2 : M * 1.1 and a = 0° 

o 


For Case 2, calculations were performed 
using only the AF algorithm. The steady results 
were obtained by cycling the step size through a 
range of values between At = 0.01 and 1.0. The 
total lift and moment were converged after 
approximately 150 steps and the potentials were 
converged to approximately 10~6.5 after 400 
steps. The resulting steady pressure 

distributions are compared with the experimental 
data in Fig. 8. For this supersonic freestream 
case, the upper and lower surface pressure 
levels are well predicted except near the 
leading edge. In general, the pressure 
distributions computed using the AF algorithm 
are in good agreement with the experimental 
data. 


Unsteady results were obtained using the AF 
algorithm for the wing oscillating with 
amplitude aj * 0.267° at a reduced frequency 
of k = 0.116. A convergence study was performed 
using 100, 200, 300, and 400 steps per cycle 
which corresponds to At = 0.2708, 0.1354, 

0,0903, and 0.0677, respectively. Unsteady 
pressure distributions for 200 and 300 steps per 
cycle are plotted in Fig. 9 along with the 
experimental data. The upper surface pressures 
are shown in Fig. 9(a) and the lower surface 
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(a) upper surface. 



x/c 

(b) lower surface. 


Fig. 6 Convergence study using AF algorithm for 
F-5 wing unsteady pressure distributions 
due to wing pitching at M « 0.9, 

Oq = 0°, = 0.109°, and k * 0.137. 
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(a) upper surface. 
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(b) lower surface. 


Fig. 7 Comparison of F-5 wing unsteady pressure 
distributions due to wing pitching at 
M * 0.9, Oq * 0°, * 0.109°, and 

k » 0.137. 


AF algorithm 
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Fig. 8 Comparison of F-5 wing steady pressure 
distributions at M * 1.1 and oq = 0°. 


pressures are shown in Fig. 9(b). Similar to 
the unsteady pressures of Case 1, approximately 
300 steps per cycle are required in Case 2 for 
convergence. Along the upper surface, the AF 
pressures generally compare well with the 
experimental data except in the real part near 
the leading edge. Along the lower surface, the 
AF results again compare well with the data and 
the calculation predicts the large change in the 

real part near 10% chord at n = 0.51, 0.72, and 
0.88. Therefore, the AF algorithm is accurate 
and efficient for application to supersonic as 
well as subsonic freestream cases. 


Concluding Remarks 

A time-accurate approximate factorization 
(AF) algorithm has been developed for the 
solution of the unsteady transonic small - 
disturbance equation. The new algorithm was 
developed as an alternative solution procedure 
to the alternating-direction implicit (ADI) 
algorithm. The AF algorithm consists of a time 
linearization procedure coupled with a Newton 
iteration technique. For unsteady flow 
calculations, the solution procedure involves 
two steps. First, a time linearization step is 
performed to determine an estimate of the 
potential field. Second, Newton iterations are 
performed to provide time accuracy. In general, 
only one or two Newton iterations are required 
to achieve acceptable convergence. For steady 
flow calculations, Newton iterations are not 
used since time accuracy is not necessary when 
marching to steady -state. 


10 




AF algorithm 
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(b) lower surface. 


Fig. 9 Comparison of F-5 wing unsteady pressure 
distributions due to wing pitching at 
M - 1.1, oq * 0°, a i * 0.267® * and 
k - 0.116. 


The AF algorithm has been used to calculate 
steady and oscillatory transonic flows at 
subsonic and supersonic freestream conditions. 
These calculations were performed for the F-5 
wing and comparisons were made with experimental 
data to assess the accuracy of the algorithm. 
The results were further compared with pressure 
distributions computed using the ADI algorithm 
for the subsonic freestream case. In general, 
the steady pressures calculated using the two 
algorithms were nearly identical and both sets 
of results compared well with the experimental 
data. The ADI results were computed using 4000 
time steps whereas the AF results were computed 
using only 400 time steps to achieve similar 
convergence. The new algorithm therefore gives 
an order of magnitude decrease in computational 
cost for steady-state calculations. Unsteady 
pressures due to harmonic wing pitching were 
very similar for the two algorithms and again 
both sets of results agreed well with the 
experimental data. The ADI results were 
computed using 2000 steps per cycle with a time 
step based on numerical stability 
considerations. The AF results, however, were 
computed using only 300 steps per cycle since 
the step size is selected for accuracy rather 
than for stability. Hence for application to 
unsteady flow problems, the AF algorithm again 
yields a significant decrease in computational 
cost. 


For a supersonic freestream case, the 
pressure distributions computed using the AF 
algorithm were in good agreement with the 
experimental pressures for both steady and 
oscillatory flows. (An ADI solution for this 
case was not obtained since the algorithm is 
unstable for swept-wing cases with supersonic 
freestream conditions.) Similar to the subsonic 
freestream case, the steady-state AF results 
were obtained using only 400 time steps and the 
oscillatory AF results were obtained using 
only 300 steps per cycle. Therefore, the AF 
algorithm is robust and efficient for 

application to steady or unsteady transonic 
flows with subsonic or supersonic freestream 
conditions. The new algorithm can provide 
accurate solutions in only several hundred time 
steps yielding a significant computational cost 
savings when compared to the ADI method. 
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